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ABSTRACT 



Disk accretion onto a weakly magnetized central object, e.g. a star, is inevitably ac- 
companied by the formation of a boundary layer near the surface, in which matter slows 
down from the highly supersonic orbital velocity of the disk to the rotational velocity of 
the star. We perform high resolution 2D hydrodynamical simulations in the equatorial 
plane of an astrophysical boundary layer with the goal of exploring the dynamics of 
non-axisymmetric structures that form there. We generically find that the supersonic 
shear in the boundary layer excites non-axisymmetric quasi-stationary acoustic modes 
that are trapped between the surface of the star and a Lindblad resonance in the disk. 
These modes rotate in a prograde fashion, are stable for hundreds of orbital periods, and 
have a pattern speed that is less than and of order the rotational velocity at the inner 
edge of the disk. The origin of these intrinsically global modes is intimately related to 
the operation of a corotation amplifier in the system. Dissipation of acoustic modes in 
weak shocks provides a universal mechanism for angular momentum and mass trans- 
port even in purely hydrodynamic (i.e. non-magnetized) boundary layers. We discuss 
the possible implications of these trapped modes for explaining the variability seen in 
accreting compact objects. 

Subject headings: accretion, accretion disks - hydrodynamics — waves - instabilities 



1. Introduction 

One of the outstanding problems in astrophysical accretion disk theory is the physics of the 
boundary layer (hereafter BL). The BL is the innermost region of the accretion disk in which the 
rotation profile of the star attaches smoothly to that of the disk ( |Lynden-Bell fc Pringlep~974 ) . In 



the BL, the rotation velocity of the disk fluid necessarily rises with increasing radius, which has 
important implications for angular momentum transport. BLs occur in a variety of systems — 
young stars, white dwarfs, neutron stars — whenever the accretion rate M is high enough for the 
disk to extend all the way down to the surface of the central object (hereafter referred to as "star" 



1 Department of Astrophysical Sciences, Princeton University, Ivy Lane, Princeton, NJ 08540; 
rrr ©astro . pr incet on . edu 



2 Sloan Fellow 



- 2 - 



for simplicity), without being disrupted by the stellar magnetic field. In such cases, up to half of 



the accretion energy is dissipated in the BL (Kluzniak 1987), which leads to intense local heating 



and gives rise to X-ray and ultraviolet emission in accreting neutron stars (Inogamov & Sunyaev 



1999) and white dwarfs ( Pringle||1977 ; Pringle &; Savonije |1979| |Narayan fc Popham||1993t |Popham| 
fc Narayan|[l995l |Piro fc Bildsten][2004l ). 



Since the rotational velocity of the gas passing through the BL changes with radius, some 
mechanism of angular momentum transport must operate there, ultimately leading to the radial 
mass transport. The nature of this mechanism is not understood at the moment. The magnetoro- 
tational instability (MRI; |Velikhov| ( |1959[ ); IChandrasekhar] fll960D ; |Balbus fc Hawley| fll991a|bD ) 
usually invoked to explain mass transport in accretion disks does not operate in the BL because 
dti 1 /dw > there [vo is the cylindrical radius). Magnetic field carried into the BL by accreting 
gas is sheared by differential rotation ( |Armitage||2002 ; Steinacker &; Papaloizou 1 |2002[ ) but whether 
this can lead to sustained mass accretion has not been conclusively demonstrated. Other transport 
mechanisms — Kelvin-Helmholtz instability (Kippenhahn & Thomas 1978), baroclinic instability 



flFujimoto|[l993l ), Tayler-Spruit dynamo ( |Tayler||1973] |Spruit||2002HPiro fc Bildsten||2004D have also 
been proposed but whether they can operate in the BL and give rise to efficient momentum trans- 
port there is far from clear. In the absence of a good understanding of the momentum and mass 
transport, even the morphology of the BL has been a matter of debate, with the proposed geome- 
tries ranging from a BL contained inside the star( |Kippenhahn fc Thomas|1 978) to a spreading layer 
extending to significant latitudes over the stellar surface ( Inogamov fc Sunyaev|[l999 , 2010). 



The BL has also been associated with variability in accreting systems that occurs on timescales 
comparable to the dynamical timescale at the surface of the star. One example of such variability 
are dwarf nova oscillations (DNOs) in accreting white dwarf systems ( Warner||2004 ), which do not 
have a definitive explanation. DNOs are characterized by an oscillation period of Pdno ^ Pk(R*), 
where Pdno is the DNO period and Pk(R*) is the Keplerian rotational period at the surface of the 



star (Patterson 1981). Since the DNO period is comparable to the Keplerian rotational period at 



the surface of the star, it is natural to associate DNOs with the boundary layer or with the region 
of the disk directly adjacent to it. 

In this work we describe a set of first principles, two-dimensional (2D) hydrodynamical sim- 
ulations of the equatorial plane of the disk-star system carried out in cylindrical geometry. The 
main result of our work is that we generically observe the excitation of a non-axisymmetric mode 
in the BL with a unique wavelength and pattern speed that is stable for many hundreds of orbital 
periods. 



The mechanism of excitation for this mode is the sonic instability (Glatzel 1988; Belyaev &; 



Rafikov 2012), which is a type of shear instability that occurs for supersonic, highly compressible 
flows. The sonic instability is closely related to the Papaloizou-Pringle instability (Papaloizou 



&; Pringle 1984; Narayan et al. 1987) and operates very differently from the more familiar KH 



instability, which occurs in the subsonic regime. 
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Dissipation of the non-axisymmetric BL mode by virtue of weak shocks results in angular 
momentum and mass transport in the BL. It has also been argued that the presence of a non- 
axisymmetric mode in the BL or on the surface of the star can modulate stellar emission, providing 
an explanation for DNOs flPapaloizou fc Pringle|[T^|Pophan^|1999t[Prro fc Bildsten||2004D . Previ- 



ously, however, there was no acknowledged mechanism for exciting non-axisymmetric modes on the 
surface of the star. Since such modes emerge naturally in our simulations, our work may help to 
understand the mechanism of DNOs, or more generically, any periodic phenomena having a period 
P~P K (IQ. 

The paper is organized as follows: In ^2] we present the model we use in our simulations, and 
in ^3] we present our results. We find that the boundary layer is unstable to sonic instabilities 
(j |3.1| ), which tend to excite a single surface mode. This mode is characterized by vortices at the 
base of the BL and by shocks launched from the top of the BL that propagate through the disk 
and are reflected back towards the BL at a Lindblad resonance. The shocks are weak with a small 
compression ratio, and we find that their properties can be easily understood in terms of the WKB 
theory for weak disturbances in the disk (j |3.2| ). We also study the stresses in the BL caused by the 
instabilities that operate there (j |3.3||3.5| ), transport of mass (j |3.6| +), how the pattern speed of the 
modes is affected by Mach number (j |3.7| ), and how the BL thickness varies as a function of both 
time and Mach number (j |3.8| ). In jjij we simulate the full range of azimuthal angle from — 2tt 
and observe long wavelength non-axisymmetric modes. 



2. Numerical Model 



We perform a set of 2D hydrodynamical simulations in the equatorial (r — <j)) plane of a star 
with an adjacent accretion disk. This setup ignores the vertical structure of the disk and star 
but allows us to capture the formation of the long-lived, non-axisymmetric structures in the BL 



at reasonable numerical cost. Similar to Armitage (2002) and Steinacker & Papaloizou (2002) we 



disregard thermodynamical evolution of the system by adopting the isothermal equation of state 
in our simulations, so the pressure is given by 



Us 2 



(1) 



Here X is the density, and s is the isothermal sound speed, which we assume to be the same through- 
out the simulation domain. Adoption of isothermal equation of state is equivalent to assuming fast 
cooling in the BL. This may not be realistic for all systems, but we still adopt this approximation 
since we are primarily interested in providing a proof of concept for the trapped acoustic modes, 
which are the main focus of our paper. 



Previously, Armitage (2002); Steinacker & Papaloizou (2002); Romanova et al. ( 2011[ ) have 
performed 3D MHD simulations of the BL. We plan to include magnetic fields in future 3D sim- 
ulations, but first it is useful to investigate and understand the restricted 2D hydro case, since it 
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allows us to study in detail trapped modes excited by sonic instabilities, which are the main focus 
of our paper. 

We numerically evolve the following set of equations 

8Y 

— + V-(St;) = > (2) 
+ V • CZvv) + VP + = 0, (3) 

at 

in the equatorial plane of the disk+star system, where $ is a fixed, time-independent potential. 
The system of equations is closed by the isothermal equation of state, equation 0. 

For the numerical evolution of equations ^ and ([3]), we use the Godunov code Athena (Stone 



et al||2008 ) in cylindrical coordinates ( Skinner fc Ostriker||2010 ). Athena is highly-parallelized and 



exhibits very low levels of numerical viscosity (j |3.6| ) making it well-suited for studying the long-term 
evolution of the boundary layer. 

We use a significantly higher resolution than previous authors, which lets us properly resolve 
the scale height of the star. This is important for studying modes excited in the BL, since the scale 



height represents a natural length scale in the problem, palsara et al. (2009) have also run high 



resolution 2D simulations of the BL in the meridional (r — 9) plane, but the axisymmetric setup 
of their simulations precluded them from seeing the sonic instability. Moreover, their simulations 
were run for ~ 1 orbital period, whereas we run our simulations for hundreds of orbital periods, as 
measured at the inner radius of the disk. This allows us to investigate the long-term evolution of 
the BL and observe the stability of modes excited on the star on very long timescales. 



2.1. Physical Setup 

In our initial setup, we consider a non-rotating star in hydrostatic equilibrium which is sur- 
rounded by a rotationally supported disk of constant density. We nondimensionalize quantities by 
setting the radius of the star to R* = 1, the Keplerian orbital velocity at the surface of the star 
to vk(R*) = 1, and the surface density in the disk to S = 1. Thus, time is scaled such that the 
Keplerian orbital period is Pk = 27r at w = 1, where w is the cylindrical radius. We take the 
potential in the system to be given by the fixed cylindrically-symmetric potential 

*(w) = -1/G7, (4) 

and we ignore self-gravity in the simulations. Since the potential given by Equation|4]is cylindrically 
symmetric, our calculations do not include the vertical stratification of the disk. 

At the start of the simulation, the star is joined to the disk via a thin region of width Sbl,o <^ 1? 
inside of which the rotational velocity rises very nearly linearly from zero (the velocity in the 
star) to the Keplerian orbital velocity (the velocity in the disk). The initial rotation profile is given 
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w < 1 - %^ "star" 



fi(ts7) = < zu~ 3 / 2 ( + i ) 1 - %^ < 07 < 1 + %^ "interface" , (5) 



by 



" 3 / 2 B7 > 1 + %^ "disk" 

and from now on, when speaking about the initial conditions of the simulations, we will refer to the 
region m < 1 — 5bl,o/2 as the "star", the region 1 — 5bl,o/2 < m < 1 + 5bl,o/2 as the "interface", 
and the region m > 1 + 5^,0/2 as the "disk", see Figure [lj It is important to bear in mind, though, 
that the rotation profile will evolve in time as a result of sonic instabilities (discussed below) that 
induce angular momentum transport. These instabilities will rapidly cause the initial "interface" 
region to thicken into a self-consistent BL, see §3.1| 



The initial density profile is specified everywhere throughout the domain through the equation 
of hydrostatic equilibrium 

1 dP _ d$ + ^ 
S dm dm 

This results in S = const within the disk, and S oc exp[— &(m)/s 2 ] within the star. Figure [l] 
shows the initial rotation profile from equation ([5]), as well as the initial density profile which is 
determined by hydrostatic equilibrium. 




Fig. 1. — Initial rotation profile Q(w) (equation [5], solid line), and the logarithm of the initial 
density (equation (6j, dashed line) for a simulation with M = 6. Note that the jump in at m = 1 
is not discontinuous and is resolved in the simulations. Various regions of the simulation domain 
are indicated. 

For each of our simulations we define a characteristic Mach number M = where s is the 
isothermal sound speed. In our units the gravitational acceleration is g = 1 at m = 1 (equation El), 
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which means that M corresponds to the actual Mach number at w — 1 for a Keplerian disk, since 
the Keplerian velocity is vk = 1 at w — 1 in our simulations. We also define a characteristic radial 
scale in the star using 

h a = M~ 2 R* = M~ 2 . (7) 
This is the exact pressure and density scale height at vo = 1 for an unrotating flow since 

(8) 



s 2 



and g can be written in terms of the Keplerian velocity as g = v 2 K jw. Plugging this into equation 
Q and using vo = 1, we recover equation Q. Note that h s is different from the scale height h d of 
the accretion disk in the vertical direction (which is not considered in our simulations). The latter 
is given by h d = s/n K so that h s = h 2 {R*)/R* = h d {R*)M- 1 < h d {R*). 

At this point it is clear that if the isothermal equation of state is assumed, the only free 
parameter in the problem is the Mach number. However, we emphasize this point by explicitly 
showing that this is the case. To do so, we non-dimensionalize the isothermal fluid dynamics 
equations by scaling the unit of length according to 1?*, the unit of time according to and 
the density according to So, which is some arbitrary density, such as the density of the disk which 
is a constant in our model. Assuming a Keplerian rotation profile {\jw form of the potential), 
equations ^ and ^ become 



<9(SV 
dt' 



^ + V ■ (SV) = 0, (9) 
+ V' • (S W) + M _2 V / S / + S^- 2 = 0, (10) 



where £ = vo / 'it!*, and the primes denote the non-dimensionalized forms of the variables and opera- 
tors, i.e. S / = S/So, V 7 = i?*V, etc. When the equations are written in non-dimensionalized form, 
it is clear that the Mach number is the only free parameter in the problem for a fixed geometry. 
Thus, the solutions for the isothermal boundary layer form a one-parameter family in the Mach 
number, and the Mach number itself plays a role which is analogous to the Reynolds number for 
the incompressible Navier-Stokes equations. Of course, if the isothermal assumption is relaxed, the 
solution will also depend on other dimensionless parameters. However, since this paper is a proof 
of concept, we do not concern ourselves with those parameters here. 



2.2. Numerical Details 



We now summarize the numerical details which pertain to all of our simulations, and we begin 
by discussing the condition of hydrostatic equilibrium. Zingale et al. (2002) discuss how to set the 
initial conditions to achieve numerical hydrostatic equilibrium in a Godunov code. However, we opt 
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for a simpler approach in which we initialize the simulation according to a numerical integration of 
equation ^ and wait for the fluctuations from the initial transient to damp out. 

The amplitude of the initial fluctuations depends on the resolution of the simulation, but 
even small fluctuations close to the inner radial boundary can amplify to become shocks as they 
propagate into the disk. This is due to the fact that the amplitude of a wave propagating in an 



isothermal atmosphere is oc p x / 2 (Landau & Lifshitz 1959) and there is a large density contrast 



between the density in the disk and the density at the inner radius of the simulation domain (factor 
of <~ 10 6 ). For this reason, we burn in all of our simulations by running the analytic hydrostatic 
equilibrium solution for a time of t ~ 100 with no perturbations, which allows the simulation to 
settle down to a numerical hydrostatic equilibrium. 
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Table 1: Simulation label, dimension of the box in cells, w & <j) extent of the box, Mach number. 

We then restart the simulation by adding a random perturbation to the radial velocity, Sv m ~ 
0.001, at each grid point for vo > 1, and for the remainder of the paper we take t = to be the 
time when the perturbations are added. The reason not to create perturbations to the velocity 
within the star stems from the fact that the density rises exponentially inside the star, and any 
disturbance propagating from the star towards the disk amplifies exponentially as well, with the 
relative amplitude going as the inverse square root of the density ( Landau fc Lifshitz||1959 ). 

The boundary conditions we use for the simulations are periodic in (j) and "do-nothing" at 
the inner and outer radial boundaries. The do-nothing boundary condition simply means that the 
fluid variables on the boundary retain their initial values for the duration of the simulation. This 
is a convenient boundary condition to use, since essentially no action takes place near the inner or 
outer boundaries: there is no mass inflow or outflow through these boundaries, and any incident 
waves are largely absorbed with minimal reflection. 



In order to accurately capture the physics of modes excited on the surface of the star (j |3.2| ), it is 
both necessary that the simulation extend to several stellar radii, and that we resolve the radial scale 
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height (equation [7]) inside the star. To satisfy the second of these conditions, we use ~ 15 — 30 cells 
per radial scale height h s within the star, depending on the simulation. However, the first condition 
that the simulation domain extend to several stellar radii makes it computationally prohibitive to 
run the simulations at realistic Mach numbers, since the scale height goes as h s = M~ 2 . As a 
compromise, we use M — 6 — 15 for our simulations, even though this is somewhat low in an 
astrophysical context. For example, a typical accreting white dwarf has an inner disk temperature 
of T - 10 5 K, a stellar radius of R* ~ 10 9 cm, and a mass of ^0.6 Mq, meaning that 

V R, kT) \0MI e R, T I x ' 

for such a system. 

In all of our simulations we use 5bl,o — 0.01 for the width of the interface, so that 5bl,o ~ h s . 
This is as thin as we can make the interface region, while still being able to resolve the physics that 
goes on there. 

Table [l] summarizes the parameters of our simulations. In all of our simulations, we only vary 
one physical parameter which is the Mach number. Simulations that have different Mach numbers 
are labeled by a different letter, whereas simulations that differ only in their numerical parameters, 
such as resolution or box size, are labeled by the same letter, but by a different number. 



Results 



3.1. Sonic Instability 



Although hydrodynamical disks with a Keplerian rotation profile are stable to axisymmetric 
perturbations by the Rayleigh criterion, the large shear initially present in the interface (equation 
[5]) makes that region unstable to non-axisymmetric shear instabilities. Since the initial flow has a 
supersonic drop in the velocity across the interface, the instability that sets in is not the classical 



KH instability for an incompressible fluid, but rather the sonic instability (Glatzel 1988; Belyaev 



& Rafikov 2012[ ). This is a global instability that cannot be derived from a local analysis and 



is similar to the Papaloizou-Pringle instability (Papaloizou h Pringle 1984; Narayan et al. 1987 



Glatzel|p88| ) 



The sonic instability operates in one of two ways, either by overreflection of sound waves from 
a critical layer (where the radial wavenumber becomes equal to zero) or by radiation of energy away 



from the BL, see |Narayan et ah] ( |1987[ ); |Glatzel| ( |1988[ ); |Belyaev fc Rafikov| ( |2012[ ). For a rotating 
flow, the critical layer of a mode is equivalent to the corotation resonance, which is located at the 
radius wqr where Vl{wcr) equals the pattern frequency, ftp. Modes, which have a corotation 
resonance (critical layer) inside the fluid domain are susceptible to instability. 



The origin of the sonic instability is closely related to the existence of a conserved action 
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(or pseudo-energy) associated with the mode, which is quadratic in perturbed fluid variables and 
changes sign at corotation ( Narayan et al.|[T987 ). A wave incident upon one side of the corotation 
region with a positive sign of the action is predominantly reflected. However, there is partial 
tunneling through the evanescent region around corotation (see j |3.2| ), and the wave emerging on 
the other side of wcr has the opposite sign of the action (i.e. negative). Global conservation of 
action then requires the reflected wave to have greater positive action than that of the incident wave 
implying amplification of the wave on that side of corotation upon reflection. If a wave is trapped 
between two corotation resonances or a corotation resonance and solid wall, this amplification 
mechanism, known as a corotation amplifier (Mark 1976; Goldreich fc Tremaine|[l978 ) , ultimately 
results in instability. 

An important difference between the sonic and KH instabilities is that the sonic instability 



operates regardless of the density contrast between the star and the disk. Belyaev & Rafikov (2012) 



have investigated the operation of the sonic instability in the presence of a density discontinuity 
for a linear shear profile. They found that the growth rate for the sonic instability was almost 
independent of the density contrast across the discontinuity. Indeed, even in the case when the 
shear flow occured over a plane reflecting wall, which can be thought of as the limit of infinite 
density contrast, the growth rate was comparable to the case of no density contrast at all. 

In our simulations, the signature of the sonic instability during the linear stage of growth is 
the presence of sound waves that are generated within the shear layer. In our initial setup, the 
region that has the greatest shear is the interface, and Fig. [2] shows the radial velocity, at 
t = 6 for simulation CI. One can compare this with Fig. 5a of |Belyaev fc Rafikov (2012), which 
shows the linear stage of the sonic instability in a uniform density medium for a shear layer with a 
piecewise linear velocity profil^j The similarities between the two figures confirm that our initial 
setup is unstable to the sonic instability. In particular, a critical layer where the pattern speed 
of the mode matches the azimuthal speed of the flow and the perturbation undergoes a phase 
shift in the azimuthal direction can be seen in Figure [2] close to the top of the interface region, at 
w « 1.005. Such a critical layer is a clear signature of sonic instabilities and is necessary for them 
to operate. We note that the sound waves emitted from the interface appear to vanish in Fig. [2] as 
they propagate inwards. This is due to the exponentially rising density towards the interior of the 
star, which causes their amplitude to decay as X -1 / 2 . On the contrary, there is no stratification in 



Fig. 5a of Belyaev & Rafikov (2012) so sound waves propagating away from the shear layer do not 



diminish as quickly. However there is still some reduction in amplitude due to the fact that sound 
waves emitted at later times have a larger amplitude than sound waves emitted at earlier times, 



due to the presence of an instability (see §4.2 of Belyaev h Rafikov (2012)) 



Belyaev & Rafikov (2012) have also shown that a high numerical resolution, of order several 
hundred cells across the shear layer for M — 10, is necessary during the initial stages of evolution 
to correctly capture the linear growth rate of the sonic instability. Because it would be prohibitive 



1 Fig. 5a of 



Belyaev & Rafikov 



(2012) is rotated by 90 degrees relative to Figure 
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Fig. 2. — a) w — <j) plot of the radial velocity at t — 6 for simulation CI (M = 12) in the vicinity 
of the interface, showing the development of sonic instability. Sound waves propagating in the flow 
both into the star (for w < 1) and into the disk (for w > 1) are clearly visible. The black arrow 
shows the direction of the unperturbed flow. 

to use such a high resolution in our simulations (we use ~ 30 cells across the interface), we do not 
attempt to accurately resolve the growth rate of the sonic instabilities during the linear regime. 
Despite this, we show in §3.8| that simulations with different numerical resolutions converge at late 
times, so the resolution does not affect the long term evolution of the BL. Moreover, using a lower 
resolution underestimates the growth rate, so going to higher resolution would simply make the 
initial instabilities proceed faster, rather than stabilizing the instability. 

After t ~ 10 the sonic instabilities saturate and after t ~ 200 the simulations settle down to 
a quasi-steady state. In this quasi-steady state, the original interface has thickened substantially 
forming a BL. For simplicity, we define the self-consistent BL to be the region in which 

0.1 < (v^(m))/v K (R*) < 0.9, (12) 

where vk(R*) — 1 is the Keplerian velocity at the surface of the star and {v^{w)) is the azimuthally- 
averaged azimuthal velocity. 

3.2. Trapped Modes 

By the time a quasi-steady state has been established and the sonic instability has saturated, 
large scale vortices are present at the base of the BL. Fig. [3] shows an image of v m at t = 298 for 
simulation A2, and the vortices are clearly visible at the base of the BL. We mention that in this 
particular case, the number of vortices is set by the box size in the azimuthal direction, but later 
in Q we discuss simulations that span the full 2tt in azimuthal angle for which this is not the case. 



radial velocity 
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Associated with the vortices is a deformation of the interface, which is shown in Fig. [3] by the 
curved streamline. Since the flow in the disk is supersonic over the surface of the BL, information 
about the interface deformation cannot propagate upstream. As a result, oblique shocks are created 
that guide the supersonic flow in the disk around the deformed surface of the BL. These shocks are 
clearly visible in Fig. [3] at radii w > 1. 

For each bump in the BL, there are two weak shocks having 5S/(S) ~ 0.1 — 0.2 < 1, where 



(E(C7)) = — J #E(C7,< 



2ir 



is the azimuthally-averaged density, and 



5E = E - (E) 



(13) 



(14) 



is the density perturbation. One of the shocks is an outgoing shock, and the other is an incoming 
shock, which is simply an outgoing shock that has been reflected from a Lindblad resonance within 
the disk (we discuss the details of this process below). The incoming and outgoing shocks intersect 
within the disk, which leads to shock crossings that result in a local enhancement of the density 
perturbation (Fig. [4^). The whole structure consisting of the vortices, deformed interface, and 
shocks rotating in a prograde fashion with a pattern speed < ftp < 1. 

We now consider in more detail the structure of intersecting shocks caused by supersonic flow 
of disk material over the surface of the star. Since the shocks are weak, we can apply linear theory 
to study their properties. In the linear approximation, the dispersion relation for a normal mode 
perturbation of the form 



is given by 



5T< oc exp 



k(w) 



k(zu')dzu' + m((p — flpt) 



-y/m 2 [to{w) 



(15) 



(16) 



which is simply the WKB dispersion relation for a fluid disk in the absence of self-gravity (Binney 



& Tremaine 2008). Here, k{w) is the radial wavenumber, rn is the azimuthal wavenumber, ftp is 



1/2 

the pattern speed, and k{vd) = (4fi 2 + wdfl 2 /dm) is the epicyclic frequency. 



For a given pattern speed ftp, the right hand side of equation (16) consists entirely of known 



quantities so the radial wavenumber as a function of radius is uniquely determined. There are two 
corotation resonances in the system, one of which is located inside the boundary layer with the 
other one located inside the disk. Each of the two corotation resonances is flanked by two Lindblad 
resonances which occur at ft = flp±K,/m. In between a pair of Lindblad resonances, k is imaginary 
(equation [l6]), and these regions are forbidden to propagating waves, i.e. the waves are evanescent 
there. 
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Fig. 3. — r — (j) plot of the radial velocity v w in the vicinity of the BL at t — 298 for simulation A2 
(M = 6). The general sense of the flow is from bottom to top, and the black lines are streamlines 
which trace the deformation of the interface. A pair of vortices are seen at the base of the BL at 
w « 0.9, and two pairs of weak shocks, each consisting of a leading and a trailing shock, are seen 
originating from w « 1. 



It should be noted that the WKB approximation, which assumes kzu/m ^> 1 breaks down 
near the Lindblad resonances, since k = at a Lindblad resonance. Nevertheless, we find that the 
WKB approximation works well in the disk, even near a Lindblad resonance. On the other hand, it 
performs poorly within the boundary layer, perhaps because the WKB approximation ignores radial 
density gradients which are large within the BL. For this reason, we will use a Keplerian rotation 
profile, = hj -3 / 2 , rather than using profiles for k and Vt determined from the simulations. A 

Keplerian profile is accurate inside the disk, and using it has the advantage of simplifying some of 
the analysis without compromising on the physics. For a Keplerian rotation profile, the corotation 
resonance in the disk is located at 

^cr = ^p 2/3 , (17) 
and the Lindblad resonances are located at 

/ \ -2/3 

zplr = ( — XT^P ) • (18) 
When a shock launched from the BL encounters the forbidden region in the disk between the 
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Fig. 4. — (a) Plot of relative density perturbation 5E/(E) at t = 284 for simulation A2. The density 
is enhanced at shock crossings, but otherwise 5E/(E) ~ 0.1 — 0.2. (b) Plot of radial velocity v w at 
t = 284 for simulation A2 (M = 6). The black curve shows the analytic solution for the shock front 
from the dispersion relation ( fl6| ) using Qp = 0.335. The dashed vertical lines show the locations 
of the inner and outer Lindblad resonances, and the solid vertical line shows the corotation radius. 
The black arrows show the direction of the unperturbed flow. 



two Lindblad resonances, it is partially reflected back towards the BL. Moreover, this reflection 
is nearly perfect, and the transmission coefficient through the forbidden region is small - so small 
that it cannot be measured in the simulations. The reflected shock then propagates inward and is 
reflected again within the BL. The reflection within the BL could either be due to the presence of 
Lindblad and corotation resonances within the BL or because the radial scale height within the BL 
becomes smaller than the radial wavelength of the shock, i.e. kh s < 1. This process of reflections 
repeats cyclically resulting in a trapped mode between the BL and the forbidden region in the disk. 

We point out that although the trapped modes are acoustic in nature, they are distinct from 
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the sonic instabilities described in §3.1| First, the sonic instabilities saturate very early in the 
simulations (t <~ 10), whereas the trapped modes only develop much later around t ~ 300. The 
likely reason for this is that the boundary layer must settle down from the violent initial state that 
occurs when the sonic instabilities saturate before the trapped modes can develop. Second, there 
is clear leakage of wave action past the critical layer during the sonic instability phase (e.g. Figure 
[2]), which is unbalanced by the wave dissipation, resulting in growth of the wave amplitude. At the 
same time, during the trapped mode phase wave action tunneling past the corotation region is not 
so noticeable (e.g. Figure [I]). As a result, a quasi-steady state is established during this phase of 
evolution, and the trapped modes are stable for hundreds of orbital periods. 

The black line in Figure [3b is the solution for an outgoing shock according to the dispersion 



relation (16) assuming a Keplerian rotation curve for ft(w) and using Qp = 0.335. Due to the 
periodic boundary conditions in the <j) direction, the shock exits from the top of the box and comes 
out again from the bottom. From Figure^, it is clear that the analytic solution given by the black 
line traces the numerical solution for an outgoing shock front very well, confirming the validity 
of the WKB approximation. As a sanity check we confirmed that the value of ftp = 0.335 was 
consistent with a movie of the simulation. 

The dashed vertical lines in Fig. show the location of the inner and outer Lindblad reso- 
nances, which are located at w — 1.99 and w — 2.16 respectively, assuming Vlp — 0.335. The solid 
vertical line is the corotation resonance, which is located at w — 2.07. In between the Lindblad 



resonances, k in equation (16) is imaginary. This region is forbidden to propagating waves, and 
outgoing waves incident upon it are reflected back towards the boundary layer. The outgoing shock 
(yellow) traced by the black line in Fig. [4)3 is indeed seen to be reflected into an incoming shock 
(blue) in the vicinity of the inner Lindblad resonance. 



3.3. Conservation of the Angular Momentum Flux 



A prediction of linear theory that we confirm in our simulations is conservation of the angular 
momentum flux. In its most general form, the angular momentum flux for a disk with no self-gravity 
is given by 

p2tt 

C L {w) = w 1 \ d^v+vn. (19) 
Jo 

Proof that Cl(hj) — const, i.e. that the angular momentum flux is conserved for stellar disks was 



first given by Toomre 


(1969) 


Bretherton & Garrett 


(1968) 



For small perturbations away from axisymmetry (linear theory) , the angular momentum flux is 
a second order quantity, and equation (19) becomes ( Binney fc Tremaine|2008 ; Balbus &; Papaloizou 
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1999[ |Steinacker fc Papaloizou||2002l ) 



C L (w) = 27tm 2 {Y i )5v^v VJ . (20) 



Here (S) is the average surface density (equation [13]), and the overbar denotes a density-averaged 
quantity 



It d ^ 

The velocity perturbations Sva, and Sv w are given by 



f(m,t)= J0 n2ir ;J . (21) 



Sv m v^-v^ (22) 
8V4 v^-v^. (23) 

We now test the linear theory during the quasi-steady state when there is a single mode 
rotating at a fixed pattern speed (£ |3.2| ). During the quasisteady state, we have outgoing shocks 
that undergo almost perfect reflection into incoming shocks at the Lindblad radius in the disk, so 
that Cl{^) ~ 0. Thus rather than measuring Cl(cu), we instead measure the quantity 



C L {w) = 2kvo 2 (Y)\5v^v vj I (24) 

where the vertical lines denote taking the absolute value. Clearly, Cl{w) > for nonzero pertur- 
bations. Moreover, away from shock crossings, when the outgoing and incoming shocks at a given 
radius are well-separated in azimuth, we have Cl(iu) ~ const. This is because for well-separated 
shocks, Cl(cu) measures the sum of the absolute values of the angular momentum flux carried by 
the individual outgoing and incoming shocks; the angular momentum flux is individually conserved 
for each of these shocks, implying that the sum of the absolute values of the angular momentum 
flux over all the shocks is also conserved. 

In Appendix [Aj we demonstrate that sufficiently far from the Lindblad resonances (when 
|f2(Tu) — Qp\ > hc) and away from shock crossings, Cl is related to the surface density perturbation 
(equation [ll]) via 

3 f27T 

c ^ = W^W)L * <<s)2 - <25) 

A similar result was previously obtained in the shearing-sheet approximation by |Goodman fc 



Rafikov| ( |200l| ). 



Figure [5] shows Cl from simulation A2 at t — 284, computed using equation (25) (solid line) 
and equation ( [24] ) (dashed line). It is not a problem that the simulation domain extends only from 
—0.4 < (j) < 0.4 rather than the full 2tt in azimuth, since the boundary conditions are periodic and 
we can stack domains azimuthally. The arrows in the figure indicate the radii at which the two 
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forms of Cl are approximately equal. We have checked using Figure |4p that these radii correspond 
to regions which are in between shock crossings. In these regions, the individual shocks are well- 



separated in azimuth, which is precisely the regime under which equations (24) and (25) must give 



the same answer. When the shocks are not well-separated, equations (24) and (25) no longer agree 



and the curves exhibit oscillations due to the presence of shock crossings. However, the two most 



important points are that away from the shock crossings: (1) equations (24) and (25) give similar 
values for Cl, and (2) Cl is roughly constant in radius with perhaps a slight, decreasing trend. 



x10" 




Fig. 5. — Plot of Cl(?u) at time t = 284 for simulation A2 (M = 6) computed using equation (25) 



(solid line) and equation (24) (dashed line). The arrows indicate the radii which are in between 



shock crossings and where there is good agreement between the two forms for Cl- 

We note that the angular momentum flux is not expected to be exactly conserved as the 
waves propagate since the shocks are not adiabatic and lose energy due to dissipation. However, 
accounting for the effect of dissipation on the angular momentum flux as a function of radius is 
made complicated by the fact that we have both outgoing and incoming (reflected) shocks in the 
simulation. 



3.4. Effective Value of a 



It is common practice in semianalytic studies of the boundary layer to parametrize the angular 



momentum transport using an ad hoc prescription. Kley &; Hensler (1987); Popham &; Narayan 



(1995); Inogamov fc Sunyaev (1999); Piro fc Bildsten (2004); Balsara et al. (2009) all used modified 



forms of the standard Shakura & Sunyaev (1973) a- viscosity prescription for disks, but the details 



of the prescription vary from one author to the next. The reason for such a lack of consensus comes 
from the fact that a simple a- viscosity prescription leads to supersonic infall velocities in the BL, 
which Pringle ( 1977[) has argued are unphysical on the basis that the disk should remain causally 



connected to the star. Consequently, it is not clear how the angular momentum transport should 
be parametrized inside the BL, or whether it can be parametrized at all with a simple prescription, 
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given the complicated shock structures that can develop. 

In the purely hydrodynamical case, the starting point of any prescription for the angular 



momentum transport should be the Reynolds stress. We define the Reynolds stress as (Balbus & 
Papaloizou||1999 ; Steinacker &; Papaloizou 2QQ2| ) 



T Re (m,t) = (^)5v w 5v^. (26) 



From this, it is possible to define a dimensionless parameter 

<E>* 



a Re (w,t) = 7 ^^. (27) 



It is clear that a Re is a dimensionless Reynolds stress, and this definition is consistent with the 



original definition of Shakura & Sunyaev (1973). However, one major point to keep in mind is 



that one should expect a Re to be negative inside the BL. This is because the rotation profile rises 
in the BL (dfl/dvj > 0) so angular momentum is transported inwards, spinning up the star. An 
inward transport of angular momentum, leads to a negative value of the Reynolds stress, and hence 
a negative value of a Re . 

Figure shows a time-radius image of a Re from simulation A2. It is evident that a Re is 
negative and is spatially nonzero only in the vicinity of the BL {w ~ 1). The latter is due to the 
fact that we have no MRI in the disk and that shear instabilities only operate in the vicinity of the 
BL. It is also clear that a Re is not constant in time, exhibiting temporal spikes at t = 150, t = 600, 
t = 1000, and t = 1400. We explain the reason for these spikes in §3.5| Typical maximum values 



of a Re as a function of radius during the low and high accretion states (see £3.5 for the definition 
of low and high accretion state) are a Re ~-2x 10 -3 and a Re ~ —0.05 respectively. 

Figures [6)3, c show time-radius images of the radial velocity, u^, and the mass accretion rate 
M (equation [33]) within the disk. The radial velocity in the vicinity of the BL is predominantly 
negative and the mass accretion rate positive, indicating accretion onto the central object, and there 
are temporal spikes in the radial velocity and mass accretion rate that are coincident with the spikes 
in a Re . In the disk proper outside the BL, the radial velocity exhibits alternating negative and 
positive regions, which are the imprints of the inward- and outward-propagating modes. In code 
units (Q, the typical values for the mean radial infall velocity during the low and high accretion 
states are ~-2x 10 -4 and ~ —0.01 respectively, and typical values for the mass accretion 
rate during the low and high accretion states are M ~ 10 -3 and M ~ .05, respectively. 

The negative value of the Reynolds stress in the BL yields a negative value of a Re . It is 
straightforward, however, to define a parameter, a^, which one might expect to be positive in both 
the BL and the disk. We can do this through a prescription for the kinematic viscosity 

v(zu,t) = Oiyslt, (28) 

where l t is the characteristic length scale for turbulent motions. The viscous stress is given by 

r(w,t) -{T,)uwdTl/dw, (29) 
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Fig. 6. — (a) Time-radius image of the dimensionless Reynold's stress parameter otR e 
equation (|27| for run A2 (M = 6). (b) The radial velocity, v^, for the same run. c) 



accretion rate, M, as defined in equation (33) for the same run. The units for b) and 
code units discussed in $2l 



defined in 
The mass 
c) are the 



and in the purely hydrodynamical case that we consider here, r — r# e , since there are no magnetic 
stresses. Thus, we can write a v as 



OL v {w,t) = -- 



(T,)sltwdQ/dw 

It is clear, then, that the condition for ol v to be positive is that 

sgn(r Re ) = -sgn(d^/diu), 
and the stress should vanish at the location where dQ/dzu = 0. 



(30) 



(31) 



For the value of l t in equation (30), we adopt the prescription of Popham &; Narayan (1995) 

— 2 91 -1/2 

n\ fdP/dwV 1 



l t (w,t) = 



(32) 
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and we define {dP/dw)/P = (d{P) / dzu) / (P) , for purposes of readability. This prescription essen- 
tially sets the turbulent scale height to be the smaller of the vertical disk scale height, s/Cl, and the 
radial pressure scale height, \P/(dP/dw)\. This generalization is necessary for the BL, since the 
disk scale height defines the characteristic length scale inside the disk, whereas the pressure scale 
height defines it inside the star. 

Figures [7] and [8] show the values of otR e (equation 27 ), ol v (equation 30 ), fi, and d^l/dw at 
t = 300 (low state) and t = 1000 (high state) respectively. In each of the figures, both oiR e and ol v 
have been averaged in time to reduce the level of fluctuations. 

It is clear that aR e is negative as discussed earlier, since the Reynolds stress is negative in the 



BL. On the other hand, a v is mostly positive in the boundary layer. However, the condition (31) 
is not exactly satisfied, leading to the large spike in Figs. [7)3 and [8)3 at vo « 1.15 and w « 1.2, 
respectively. Moreover, the rotation profile in the high state at t — 1000 is very flat between 



•97 < vj < 1.07 (Fig. ^p) and even contains a change in sign of dVtjdw (Fig. pll). This again leads 



to a violation of condition (31) and leads to the densely packed spikes between .97 < zu < 1.07 in 
Fig. [8J3. 

If the non-trivial nature of the rotation profile (Figs. ^ and that we observe in our 
simplified simulations holds under astrophysical conditions, it seems unlikely that the dynamics of 
the BL can be captured with a simple prescription for a y . On the other hand, the dimensionless 
Reynolds stress aR e is quite smooth and well-behaved inside the BL (Figs. [7^ and [8^1). Thus, aR e 
could be the preferred parameter for constructing semi- analytical models of the BL. 



3.5. High and Low Accretion States 

From Figure [6] it is clear that our simulations are characterized by states of high and low 
accretion rate. The high accretion states are centered ont= 150, t = 600, t = 1000, and t = 1400 
and correspond to large absolute values of the Reynolds stress and of the radial infall velocity (Fig. 
[6]). The low accretion state is relatively calm and is characterized by a single dominant mode that 



rotates at a constant pattern speed as discussed in §3.2| On the other hand, the oscillations of the 
interface during the high accretion state are more violent, and it becomes more difficult to determine 
a single dominant mode or a unique pattern speed. Nevertheless, the shock structure seen very 
distinctly in the low accretion state is still present to some extent even in the high accretion state. 
Fig. [9] shows a plot of the radial velocity in the high accretion state at t = 1000 for simulation A2, 
and can be compared with Fig. [4)3, which shows the radial velocity during the low accretion state. 
It is clear that the high accretion state is more violent and chaotic than the low accretion state, 
but a dominant mode is still present. 

We can understand the reason for transitioning from the low accretion state to the high ac- 



cretion state by considering as a function of time. Fig. 10 shows at t = 800, 900, 1000 



and 1100. These times span the duration of the third high accretion state in Fig. |6j We see that 
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between t = 800 and t = 900, a bump develops in the azimuthal velocity profile. This leads to a 
rapid rearrangement of the velocity profile between t = 800 and t = 900, which erases the bump, 
and by t = 1100 the rearrangement of the velocity profile is complete. 

The cause of the bump that appears in the profile are likely to be the oblique shocks present 
at the top of the BL. Each time a fluid element at the top of the BL passes through a shock, it is 
slowed down, and the collective effect of these shock passages creates the bump. Once the bump has 
been created, the mechanism for erasing it is the KH instability. KH instabilities can operate near 
inflection points in the velocity profile. At an inflection point, the radial derivative of the vorticity 
is constant, making it permissible to mix neighboring fluid elements and extract energy from the 
shear in the flow. In the vicinity of the bump, the vorticity profile is very flat, which allows the KH 



instability to operate there. Fig. 11 shows a snapshot of the vorticity (V x v) in the vicinity of the 
BL at t = 1000, during the height of the high-accretion state. One can see that around w ~ 1, the 
vorticity is concentrated into structures that resemble the cats eye pattern associated with classical 
KH instabilities. 
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Fig. 9. — Radial velocity during the high accretion state at t — 1000 for simulation A2 (M = 6). 
Compare this with a snapshot of the system in the low accretion state in Figure [3| 
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Fig. 10. — Profiles of mean azimuthal velocity at t = 800 (black), £ = 900 (red), t = 1000 (blue) 
and t = 1100 (green) for run A2. 
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Fig. 11. — Image of the vorticity, V x v, in high accretion state at t — 1000 for simulation A2 
(M = 6). The cat's eyes are the red oval structures with wispy arms around w — 1 and are 
indicative of KH instability. 



A curious point is that after a certain amount of time the system stops switching between 
high and low states and stays in the low state indefinitely. For instance, in Figure [6] it is clear that 
after four high-accretion states during which the system is in the high accretion state it stays in 
the low accretion state, until the end of the simulation. The cause of this behavior is unclear, but 
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one should keep in mind that the inner part of the disk becomes more and more depleted of mass 



throughout the course of the simulation, see £3.6, If the material in the inner disk were replenished, 



as it would be in a real disk, the high accretion state might not shut off and continue to operate 
periodically. 

It may be tempting to associate the high and low accretion states with the outburst and 
quiescent states of a CV. However, such a comparison is superficial given the fact that we are using 
an isothermal equation of state and have ignored magnetic fields. Nevertheless, the modification of 
the rotation profile by shocks and the subsequent onset of the KH instability is a very interesting 
phenomenon in its own right, even if it bears no relation to the mechanism causing outbursts in 
CVs. 



3.6. Mass Accretion Due to Shocks 

We confirm that even though we have purely hydrodynamical simulations, there is mass accre- 



tion onto the star. Figure 12 shows the density at t — and at t — 2000 for runs A5 and CI. From 



the figure, it is clear that the inner part of the disk has been evacuated at t — 2000 with material 



having been accreted onto the star. In Figure 12 a,, the innermost Lindblad radius is at w = 2.03. 



One can see that the depletion of mass occurs primarily inside of this radius, which is the region 
in which the trapped modes are present. There is some depletion even beyond this radius which is 
possibly due to the generation of a pressure gradient as mass in the inner region accretes. The gain 
in the radius of the star is more prominent for run A5 as opposed to run Cl, for two reasons. First, 
the Mach number is two times lower in run A5 as compared to run Cl (6 vs 12), which means the 
radial scale height is four times smaller and the compression of accreted material larger in run Cl. 
Second, more material has accreted at t — 2000 for run A5. 

Most of the mass accretion in the simulations occurs during the high state, which is also when 
the Reynolds stress is largest. However, during the low accretion state, it is possible to "predict" 
the mass accretion rate in the disk based on the value of 5S/(S), and compare it to the "observed" 
accretion rate: 

M = -2ttw{E)v^, (33) 

where the overbar represents an azimuthal density- weighted average. 

To calculate the predicted mass accretion rate, we note that during the low accretion state 
there is a well-defined structure of shocks that rotates at a constant pattern speed. Given 5E/(S), 
one can use the fact that these shocks carry angular momentum, and because of dissipation, the 
angular momentum in the shock is transferred to the bulk flow driving accretion. In Appendix [Bj 
we derive an analytical expression for the mass accretion rate M resulting from dissipation of weak 



shocks, given by equation (B7). We find that the mass accretion rates measured in our simulations 



using equation (33) agree to order of magnitude with the predicted mass accretion rates using 
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Fig. 12. — (a) Azimuthally-averaged density as a function of w at t = (dashed line) and t = 2000 
(solid line) for run A5 (M = 6). Evacuation of the part of the disk interior to the inner Lindblad 
radius (at w — 2.03) is clearly visible, with mass ending up on the star, (b) Same as (a) but for 
run CI [M = 12). 



equation (B7). For simulation A2, the mass accretion rate is M ~ 10 -3 (in our dimensionless 
units), during the low state around t ~ 300. This explains why the disk can persist for hundreds 
of orbital periods in the low state without being depleted. We note that during the high accretion 
state, the accretion rate can be much higher reaching values as high as M ~ 10 _1 — 10 -2 in the 
boundary layer, where shear instabilities operate. 

We also point out that accretion due to numerical viscosity is negligible in the simulations. For 
instance, in simulation A2, numerical viscosity gives rise to a numerical accretion rate of M ~ 10 -6 . 
This estimate was made by running a simulation with the same parameters as simulation A2, but 
without perturbations. In the absence of perturbations, the sonic instability does not set in, and we 
maintain an azimuthally-symmetric flow for the duration of the simulation, allowing us to measure 
the numerical mass accretion rate directly. We point out, however, that such an azimuthally- 
symmetric flow is grid-aligned in cylindrical coordinates. Thus, our measurement of the numerical 
mass accretion rate likely underestimates its true value in the science runs, where radial motions are 
present. Nevertheless, the fact that the accretion rate due to numerical viscosity in the grid- aligned 
case is three orders of magnitude lower than the measured accretion rate in the low-accretion 
state of simulation A2 provides convincing evidence that numerical viscosity is negligible in the 
simulations. 



3.7. Stability of the Pattern Speed 



Given that the system can transition between high and low states, one might wonder about 
the long-term stability of the modes discussed in §3.2| We find that the pattern speed in a given 
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simulation stays constant to within a few percent over the course of hundreds of orbital periods. 
This is true even if the time interval in question contains transitions between high and low states. 
Figure [13] shows an image of the radial velocity for simulation A2 at time t = 1284. This image 
can be compared with Figure^ which is from the same simulation but was taken at t = 284. The 
pattern speed at t = 1284 is approximately ftp = 0.325, which is close to the pattern speed at 
t = 284 of ftp = 0.335. Note that the interval of time between t = 284 and t = 1284 contains two 
transitions from high accretion state to low accretion state. 



radial velocity 




radius 

Fig. 13. — Density plot of v w at t = 1284 for simulation A2 (M = 6). The black curve shows 
the analytic solution for the shock front from the dispersion relation ( fl6| ) using ftp = 0.325. The 
dashed vertical lines show the locations of the inner and outer Lindblad resonances, and the solid 
vertical line shows the location of the corotation radius. 



3.8. Width of the BL 



We plot the thickness of the BL defined by equation (12) as a function of time in Figure 14 for 
runs Al, A2, A3, Bl, CI, and Dl. The thickness of the BL greatly increases in the beginning of 
the simulation due to the sonic instabilities. It then stays approximately constant during the low 
accretion states and undergoes further changes during the high accretion states, eventually settling 
down to an approximately constant value. The fact that runs Al, A2, and A3 yield almost the 
same BL width at late times suggests that the simulations are converged, since these three runs 
differ only in their resolution. 

It follows from Figure [10] that after significant evolution has taken place — typically after 
hundreds of orbital periods — the region in which azimuthal velocity transitions from the Keplerian 
profile tov ( j ) = covers the radial interval 0.9 < w < 1.2. Thus, it extends both into the star and 
into the disk. This may have interesting implications for the morphology of the boundary layer in 
astrophysical systems. 



It is clear from Figure 14 that the boundary layer becomes thinner with increasing Mach 
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Fig. 14. — Width of the BL as defined in equation (14) as a function of time. Runs Al, A2, and 



A3 all have M = 6, run Bl has M = 9, run CI has M 12, and run Dl has M 15. 



number. In Figure [15j we plot the width of the boundary layer at t = 1500 as a function of the 
Mach number on a log-log plot. The data point for M = 6 is a simple average of simulations Al, 
A2, A3 and the points for M = 9, M = 12, and M = 15 are from simulations Bl, CI, and Dl 
respectively. The dashed line shows the best least squares fit to the data, and from its slope we 
deduce that Hbl M -L9 , where Hbl is the width of the BL at t = 1500. Since the scale height 
goes as h s = M~ 2 (equation (7]), this indicates that the deceleration of the flow in the BL occurs 
over an approximately fixed number of scale heights in our simulations (~ 7 — 8), regardless of the 
Mach number. 



4. Full Disk Simulations 

Up until now, the simulations we have discussed have had a <j) extent of 0.8 radians. In this 
section, we consider global disk simulations, where the azimuthal angle runs from to 2tt. This is 
important, since it is possible that the azimuthal wavenumber of the dominant global mode is of 
order unity. Lower azimuthal wavenumbers are inherently more interesting for explaining physical 
phenomena having frequencies less than or of order the orbital period at the surface of the star, 
such as DNOs. This is because the frequency of the mode is given by uj — mflp, and only if m is low 
enough will uj also be low enough to provide a plausible explanation for such phenomena. Thus, we 
determine the azimuthal wavenumbers and pattern speeds of the dominant modes obtained from 
simulations in which <j) goes from to 2tv. 



Panels (a)-(d) of Figure 



show images of the radial velocity from simulations A4, B2, C2, 
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Fig. 15. — Boundary layer width at t — 1500 as a function of the Mach number in our simulations. 
The dashed line indicates the best fit line through the points and has a slope of —1.9. This is close 
-2, which is the expected value of the power law index for the scaling of the boundary layer width 
with the Mach number. 



and D2, which have Mach numbers of M — 6, 9, 12, and 15 respectively. These snapshots reveal 
a rather interesting structure of fluid perturbations, which is best seen at low values of M. In 
particular, in panel (a), one can clearly discern a large-scale m — 3 global mode along with some 
high frequency features, which are due to vortices at the base of the BL. There are more than 20 
individual vortices that reside in the BL (where ft goes down from its Keplerian value in the disk 



to in the star). Each of the vortices in the BL launches a shock into the disk, as discussed in {3.2 
but now the individual shocks are superimposed on top of a global m = 3 mode. 

Thus, close inspection reveals two types of periodicities in our simulations — a large scale 
mode with relatively low wavenumber m that extends out from the BL into the disk, and a set of 
periodic vortices contained in the BL with higheij^] wavenumber m v (m v > 20 in Figure 16 a). The 



vortices at the base of the BL are only visible in panels (a) and to some extent b of Figure (16{ 
However, they are present as well in panels c and d and become apparent if the color scheme is 



saturated. The origin of the global mode is traced to a geometric resonance condition in £4.1, but 
we leave the detailed investigation of the BL vortices to future study. 



2 Note that in simulations presented in { 3 



which covered only a limited interval in <j> it is generically found that 



(and typically equal to 2) since the azimuthal extent of the domain was setting the wavenumber of the mode. 
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The pattern speed of the m = 3 global mode in panel a of Figure fl6] is Qp = 0.28. This is 
comparable to the pattern speed of the shock structure for simulation A2 discussed in £3^2, which 
has the same Mach number as simulation A4. The frequency of the m = 3 mode is mVtp — 0.84, 
and thus is lower but comparable to the orbital frequency at the inner radius of the disk. 




(c) (d) 

Fig. 16. — Images of the radial velocity for simulations A4 (panel a), B2 (panel b), C2 (panel c), 
and D2 (panel d), which have M — 6, 9, 12, and 15 respectively. The view has been zoomed in, in 
panels (b), (c), and (d) to give a better view of the modal structure in the vicinity of the BL. The 
m-numbers of the dominant global modes for panels (a)-(d) are m = 3, 13, 16, and 15 respectively. 
Solid and dashed circles show the locations of the corotation and the inner and outer Lindblad 
resonances in the disk. 
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4.1. Relation of Pattern Speed to Azimuthal Wavenumber 

As the Mach number is increased, there is a general trend for the azimuthal wavenumber of the 
large scale pattern, m, to increase as well. Figure [T7^i shows the wavenumber of the dominant global 
mode as a function of Mach number. We also plot the pattern speed as a function of Mach number 
in Figure [17) 3. There is again a general trend for the pattern speed to increase with increasing 
Mach number. 




9 12 
Mach number 



0.8 



0.6 



0.4 



0.2 



9 12 
Mach number 



15 



(b) 



Fig. 17. — (a) Plot of the azimuthal wavenumber of the global mode as a function of Mach number, 
(b) Plot of the pattern speed as a function of Mach number for the global mode. 



The dependence of the pattern speed on the azimuthal wavenumber of the global mode can 
be understood if we consider the total azimuthal angle traversed by a shock in traveling from the 
boundary layer to the inner Lindblad resonance in the disk and back. This angle is given by 



Ac/) = - J miLR dw^j - n P ) 2 - 



and can be derived by using equation (16) together with the relation 



d(f) k{w) 
dw 



m 



(34) 



(35) 



which describes the trajectory along which the phase in equation (15) is constant. Here, vjbl and 



wilr denote the boundary layer radius and the inner Lindblad radius in the disk, respectively. The 
factor of two out front comes from the fact that we are considering the total change in angle for a 
shock that starts from the boundary layer, is reflected at the Lindblad radius, and comes back to 
the boundary layer. 



Formally we need to set vjbl equal to the location of the Lindblad resonance inside the BL. In 
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practice, though, the shocks in our simulations are always observed to originate very near w = 1, 
and thus we simply set wbl — 1 for the purposes of our calculations below. 

An interesting question to ask is, given that we observe a mode with azimuthal wavenumber 



m, what pattern speed would we predict based on equation (34), and how well does this predicted 



pattern speed match the observed pattern speed in our simulations? Assuming that the change in 



azimuthal angle for one "out and back" trip for the shock is given by equation (34), the condition 
that the shock pattern be closed is given by 

A*=^, (36) 

where p and q are positive integers. It is natural to assume that p = 1, especially since the global 



modes don't appear to overlap in Fig. 16 as they would if p > 1. For g, an obvious choice to make 



is q = m, since this naturally gives an m-fold pattern. Thus, we postulate that 

A0=-, (37) 
m 

where m is the azimuthal wavenumber of the global mode. 

For simplicity, we assume a Keplerian rotation profile, in which case our estimate for Acf) 
becomes 



O r^ILR / / 777-3/2 N 

A(p = - I dw^ (w-z/ 2 - n P y - [ ) , (38) 



where zujlr is given by equation (18). Assuming a Keplerian rotation profile gives a value of Acf) 
which is slightly larger than observed in the simulations. For instance, see Fig. and note in 
particular that around w ~ 1, the shock profile becomes less steep than the analytic prediction 
using a Keplerian profile (black curve), due to the fact that the angular velocity is sub-Keplerian 
in that region. However, using the values of Vt{w) and k(vj) is complicated by the fact that the 
WKB approximation breaks down in the boundary layer, (see j|3.2|for a more detailed discussion). 



The open circles in Fig. 18 show the value of f2p that gives the correct value of Acf) (equation 
37) using the formula (38) for the m- numbers observed in the simulations. Conversely, the solid 
circles in Fig. 18 show the values of ftp as measured from the simulations. We find that the 
predicted value of tip is within ~ 10% of the value observed in simulations and is always larger 
than the observed value. The latter point can be explained by the fact that using a Keplerian 
profile results in a larger value of Act> as explained above. A smaller value of Acfi would lower the 
predicted value of fip, bringing the analytical predictions into closer agreement with the simulation 
results. 

The good agreement between the predicted and observed pattern speeds over a range of Mach 
numbers suggests that there is a direct relation between the m-number of the global modes in the 
disk and the pattern speed. This relation is essentially the combination of a periodicity condition 
and the dispersion relation for the modes. 



1 
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Fig. 18. — Plot of the pattern speed of the global mode as a function of the azimuthal wavenumber 
m. The closed circles show the pattern speed measured from the simulations and the open circles 
show the predicted pattern speed (see text after equation ([38])) . 



5. Discussion 

The primary result of this work is the identification of a novel mechanism of angular momentum 
and mass transport in the boundary layer relying on non-axisymmetric acoustic modes which are 
generically excited in highly supersonic shear layers. Dissipation of these modes in weak shocks 
drives the evolution of the system and ultimately results in the mass transfer from the disk into the 
star. Our work is exploratory in nature and necessarily makes a number of simplifying assumptions, 
many of which are discussed in §5.2| Nevertheless we expect the main qualitative features of the new 
transport mechanism to persist also in more realistic settings typical for astrophysical boundary 
layers. The basic reason for this expectation lies in the fact that essentially the only ingredient 
necessary for the generation of acoustic modes is the presence of supersonic shear which is natural 
in astrophysical BLs. 

Previously a number of hydrodynamical and MHD instabilities — Kelvin-Helmholtz, baro- 



clinic, Tayler-Spruit dynamo ( |Kippenhahn fc Thomas||1978[ |Fujimoto||1993[ |Tayler||1973[ |Spruit 



2002; Piro &; Bildsten 2004; Inogamov Sz Sunyaev 2010) have been invoked to drive the angular 



momentum transport in the BL. However, none of these instabilities were actually shown to operate 
under the conditions typical for the BL, where the supersonic nature of the flow is crucial for the 
dynamics and can drastically change the way in which these classical instabilities (usually explored 
in incompressible limit) operate. Armitage (2002) and Steinacker fc Papaloizouj ( |2002[ ) have run 
direct 3D MHD simulations of BLs at a resolution that although state of the art at the time is low 
compared to the current simulations - their typical resolution in the radial direction is about an 
order of magnitude lower than in our work. These authors found that magnetic field amplification 
and mass accretion took place in the BL. However, whether the latter is due to the magnetic stresses 
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and not due to the hydrodynamical effects like the one discussed here is not obvious. Transport 
induced by acoustic modes does not require magnetic field to be present in the first place, which 
makes this mechanism quite universal. 



As part of our work we also looked at the linear stage of the sonic instability (see £3.1), 
which was previously explored under a more limited set of assumptions by Belyae v fc Ra fikov 
(2012). In particular, |Belyaev fc Rafikov (2012) neglected the Coriolis force in their calculations 
and worked in Cartesian geometry. Despite that, their main results are in agreement with our more 
general calculations run in cylindrical geometry which fully account for the effects of rotation. This 



concordance justifies the simplifications made by |Belyaev fc Rafikov (2012) in studying the sonic 
instability. 



5.1. Variability in the BL 



Because of the periodic nature of the trapped modes and their stability over many orbital pe- 
riods, one may wonder whether they can produce a periodic modulation of the disk-star luminosity. 
One natural application of such a modulation is to the variability observed during the dwarf novae 
outbursts. The leading explanation for this variability involves magnetically-channeled accretion, 



which spins up an equatorial belt on the primary (Paczyhski 1978; Warner & Woudt 2002). An 



alternate explanation, and one that is relevant to our results is in terms of modes excited on the 
surface of the WD. 



Papaloizou & Pringle (1978) found that global non-axisymmetric modes excited in the surface 



layers would have the correct periods to account for DNOs. Later Popham (1999) studied the 



possibility that DNOs could be caused by a bulge in the BL and Piro & Bildsten (2004) considered 
shallow water modes in a BL that had spread meridionally to high latitudes. However, one of the 
main impediments to such theories is identifying the physical mechanism that would excite such 
surface modes in the first place. 

If the picture of sonic instabilities leading to trapped modes that we have presented here 
remains valid even when additional physics is included (e.g. 3D nature of the flow, realistic cooling, 
effects of magnetic fields), then it would provide a mechanism for exciting surface modes. Since the 
frequency of DNOs is typically less than the orbital frequency at the surface of the star, the most 
favorable scenario for this mechanism to explain DNOs would involve a low m number azimuthal 
mode m < 3. This is because one might expect a high m- number mode to produce a small signal due 
to integrating the light over a hemisphere of the star. Also, the criterion that the DNO frequency 
be less than the orbital frequency at the surface of the star means that m cannot be too large 
otherwise the criterion cjdno — m^p < 1 will not be fulfilled. 

The only one of our simulations to satisfy the constraint mflp < 1 is simulation A4 (M = 6) 
with the higher Mach number simulations — more typical for realistic astrophysical systems — 
producing values of mVtp that are too high for DNOs. However, a detailed comparison to DNOs 
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is premature at this point given the exploratory nature of the simulations - 2D with isothermal 
equation of state. 



5.2. Extensions 

We have performed 2D hydrodynamical simulations of the boundary layer and have found ex- 
citation of non-axisymmetric surface modes by the sonic instability. One may wonder how relevant 
our results are to real 3D disks with magnetic fields. We plan to address this question with future 
simulations, but for now we comment on what we expect may change in going from 2D hydro to 
3D MHD simulations. 

We have carried out some preliminary calculations to demonstrate the fundamental difference 



between the 2D and 3D cases and show the results in Figure 19 Panel (a) of Figure 19 corresponds 
to simulation Al at time t = 1000. Panel (b) is for a 3D simulation in cylindrical geometry with 
identical parameters except: the simulation domain has a vertical extent of —0.167 < z < 0.167 
(i.e. 2hdi where hd = s/Q is the disk scale height), and there are 64 cells in the z-direction. We use 
a cylindrical potential of the form &(w) — —1/zu, where w is again the cylindrical radius, so the 
disk is unstratified. The boundary condition in the z-direction is periodic for the 3D simulation. 



The shock structure and the vortices at the base of the BL are clearly visible in Figure [19^ , 
but are absent or only barely present in Figure [Hfo . Instead, in 3D the vortices at the base of 
the BL have been replaced by turbulence. We will explore the appearance of vortices at the base 
of the BL in a future work, but we remark here that the appearance of vortices could depend on 
the aspect ratio of the BL. If the BL has a vertical extent which is small compared to its radial 
extent, then one might expect vortices to be present, whereas if the BL has a vertical extent which 
is thick compared to its radial extent, then it may instead be turbulent. Nevertheless, the structure 
of intersecting shocks is still present in the disk even in the 3D case with turbulence (Figure [TffiQ, 
although it is less clear cut than in the 2D case with vortices. One may also expect that the large- 
scale global modes discussed in ^4] would not be affected by the turbulence, since the periodicity 



of global modes is determined by the geometric resonance condition (37). This does not involve 
the physics of what goes on inside the BL, which occurs on scales that are small compared to the 
wavelength of the global modes. 

Additional effects may emerge in 3D hydrodynamical calculations in which vertical stratifica- 



tion in the disk is properly accounted for. Lubow & Ogilvie (1998) and Bate et al (2002) have 



demonstrated that waves in thermally stratified disks tend to propagate in such a way that their 
action density is concentrated predominantly near the disk surface, which increases wave amplitude 
and leads to faster nonlinear damping. In our case this effect may influence the damping of the 
trapped modes leading to their faster dissipation. Another important aspect of stratified 3D sim- 
ulations is that they would allow one to naturally explore the spreading of material on the stellar 



surface towards high altitudes, resulting in formation of a spreading layer (Inogamov &; Sunyaev 
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1999, 2010). 
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Fig. 19. — (a) Image of v m at t = 1000 for simulation A2 (M = 6). Vortices are clearly visible at 
the base of the BL. (b) Same as panel (a) but for the 3D simulation at t = 1000. The base of the 
BL is turbulent and no clear vortices are present. Nevertheless, the shock structure is still present 
in the disk, although it is not as crisp. 

One may also wonder how the addition of a magnetic field would modify our results. The BL 
is linearly stable to the MRI, since it has a rising rotation profile (dfi 2 jdw > 0). Nevertheless, 
due to the large shear present in the BL, a radial magnetic field can be kinematically sheared to 



produce a strong toroidal field. Armitage (2002) and Steinacker & Papaloizou (2002) found such 



an amplification in their simulations, but the field remained subthermal, and the magnetic energy 
density was only a fraction (~ 15%) of the thermal energy density in the boundary layer. This 
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means that even in the MHD case, most of the pressure still comes from the gas, so one might 
expect the magnetic field and its associated pressure to only provide a small correction to the 
dispersion relation of the trapped modes. 

Another way that MHD could affect the trapped modes is that the velocity fluctuations in 
the disk due to the MRI could be much larger than the velocity perturbations due to the shocks. 
Hawley et al. | (fl995) found that the velocity fluctuations due to the MRI were ~ 10% of the sound 



speed in the disk. This is comparable to the velocity perturbation induced by the trapped modes 
in the disk. Thus one might expect that the trapped modes will have their wavefronts deformed by 
the MRI turbulence, but that the overall shock structure will still survive. 



6. Conclusion 

We have performed 2D hydrodynamic simulations of the boundary layer that self-consistently 
forms when the accretion disk is brought into contact with the stellar surface. Our calculations are 
restricted to the equatorial plane of the star+disk system which does not allow us to resolve the 
disk vertically but, quite importantly, captures the non-axisymmetric structures that emerge as a 
result of evolution. In our simulations, we used an isothermal equation of state and a Keplerian 
rotation profile inside the disk. We varied the Mach number at the inner edge of the disk to study 
how it affected the outcome of the simulations. 

We find that sonic instabilities set in very rapidly in the interfacial region of large, supersonic 
shear between the star and the disk and saturate within a few orbital periods. After saturation, the 
initially thin boundary layer region expands radially both into the star and into the disk, with its 
final thickness depending on the Mach number as oc M~ 19 . This scaling implies that the boundary 
layer spans the same number of pressure scale heights (~ 7 — 8) in our isothermal simulations 
independent of the Mach number. 

After t rsj 50 orbital periods, a trapped mode develops between the BL and a forbidden region in 
the disk, which consists of a corotation resonance flanked by two Lindblad resonances. The trapped 
mode is actually a weak shock (relative density perturbation ST,/ (S) < 1), and one can think of it as 
a global acoustic wave that alternately reflects off the BL and the forbidden region in the disk. The 
origin of this mode is intimately related to wave tunneling through the forbidden region and the 
associated phenomenon of overreflection. This trapped mode can persist for hundreds of orbital 
periods at a well defined pattern speed, ftp < Qk(R*)- Both the pattern speed and azimuthal 
wavenumber of the the mode increase with increasing Mach number, and one can express the 
pattern speed in terms of the azimuthal wavenumber using a geometric resonance condition. 

The trapped modes and initial sonic instabilities are able to drive accretion onto the star via 
dissipation in weak shocks even in the absence of magnetic fields. The steady state consisting 
of a trapped mode rotating at a given pattern speed can be interrupted by high-accretion states 
that are possibly caused by KH instability in the BL due to modification of the rotation profile by 
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shocks. The non-axisymmetric nature of the acoustic modes found in our calculations gives rise to 
variability of the emission produced in the BL, which may be relevant for explaining phenomena 
such as dwarf novae oscillations in cataclysmic variables. 

We thank Jeremy Goodman for useful discussions and the referee for his comments and sug- 
gestions. The financial support for this work is provided by the Sloan Foundation and NASA grant 
NNX08AH87G. 
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A. Angular Momentum Flux from a Superposition of Modes 

In the WKB limit (k ^> m/R), the angular momentum flux for a normal mode with azimuthal 
wavenumber m becomes 



CL,mM TTSgn(k) 



mws 2 8Y^ n {vj) 



(Al) 



(see Goldreich fc Tremaine| ( |1979 ) and Appendix J of |Binney fc Tremaine (2008)). Here k is the 
radial wavenumber and sgn(fc) determines the direction of angular momentum transport, inward 
or outward. The Fourier coefficient 5S m in equation (Al) is defined such that an arbitrary density 
perturbation, 0, £), can be expressed as a sum over Fourier components having m > 0: 



8Yj(m,(f),t) — Re 8T lm (w)e 



i(mcj)—ujt) 



m=0 



(A2) 



This preprint was prepared with the A AS IAT^X macros v5.2. 
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Since the shocks in our simulations are not perfect normal modes, we need to compute the 
angular momentum flux for an arbitrary density profile with a given pattern speed. This can be 
done by summing over the angular momentum transported by the individual modes so we have 



(A3) 



m=l 



Note that the m = term does not contribute to angular momentum transport, which can be seen 



from equation (Al). This point allows us to perform a mathematical trick that will simplify the 



analysis. Since the m = mode does not contribute to the total angular momentum transported 



flux, we are free to set it to zero in equation (A2). Although this has an effect on the density 
profile, it does not change the total angular momentum flux. 



Setting the m = component of the density profile to zero, we can rewrite equation (A2) as 



oo 

E 

m=— oo 



5Y im (w)e 



im(cf)—Qpt) 



(A4) 



where T,- m {w) = S^(iu). Parseval's theorem, which we shall shortly invoke, can then be written 

as 



r 27T 



^ i pZ 

V \SE m \ 2 = - / d<j>8T? 
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(A5) 



m=l 



Using Parseval's theorem, we now derive an approximate formula for computing the angular 
momentum flux that does not involve summing over the individual Fourier coefficients as in equation 



(A3). We begin by making the assumption ft — £lp ^> n/m (k, = ft for Keplerian profile) so that 
the expression for \k\ (Eq. [l6|) simplifies to 



777 

\k\ w — (Q-flp) 
s 



(A6) 



where ft > ftp. Consequently, this assumption is true if m ^> 1, and we are far away from any 
resonances so that ft — ftp ~ ft ^ n/m. 



Substituting equation (A6) into equation (Al), we have for the angular momentum flux due 
to a given mode 



C L ,m(™) -^sgn(fc)- 



(Z)(ft-ft P )' 

The total angular momentum flux from equation ( |A3| ) is given by 

c .3 



(A7) 



Cl{™) = vr 
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(A8) 



171=1 
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We shall assume that for a given shock, all of the modes have the same sign of sgn(fe), i.e. the shock 
is either incoming or outgoing. This is a reasonable assumption except near shock crossings, where 
incoming and outgoing shocks having opposite signs for sgn(fc) interact. For a given shock, this 
allows us to take sgn(fc) out of the summation sign, assuming the individual shocks are azimuthally 
well-separated. Since Cl is the sum of the absolute values of the angular momentum flux for the 
individual shocks, assuming the shocks are well-separated, we arrive at 



3 00 

S < ro) = *<EHra < A9) 

x ' v 7 m=l 



Here we have assumed ft > ftp which is true for the trapping region between the boundary layer 



and the Lindblad resonance in the disk. Using Parseval's theorem, equation (A9) becomes equation 

ph. 



B. Mass Accretion Rate Due to Weak Shocks. 



The mass accretion rate through the disk can be expressed through the divergence of the 
angular momentum flux F w ^ s deposited in the disk fluid by the waves (or acoustic modes) due to 
their dissipation at the shock fronts ( Lynden-Bell fc Pringle|[T974 ) 



M{r) 



dl y 1 dF wAis 
dw J dr 



(Bl) 



where / = ft{w)w 2 is the specific angular momentum for circular orbits. The subsequent derivation 



is based in part on the work by Larson (1990). 



The modes are stationary in a frame rotating with angular speed ftp. Excitation of a wavepacket 
carrying angular momentum A J requires energy AE = ftpAJ. Dissipation of this wavepacket at 
some other radial distance w transfers angular momentum A J to the disk fluid — a process which 
is necessarily accompanied by the mechanical work ft{vo)AJ being done on the disk. The rest of 
energy, [ftp — fl(w)] A J, is available for heating the disk and is ultimately radiated away. 

Thus, if the dissipation of the wave adds angular momentum to the disk fluid at the rate 
dF w ^[ s /dvj per unit of w (the radial divergence of the wave angular momentum flux due to the 
wave damping), then the local tidal energy dissipation rate dE w /dw is given by 



dE w 
dvo 



[Q P -Q(w)] 



OR 



w,dis 



dw 



(B2) 



See Goodman h Raflkov (2001) for the discussion of this issue. 



It is difficult to calculate dF w ^ s /dvu from first principles. However, one can compute dE w /dva 



and use equation (B2) to obtain dF w ^ s /dzu. The reason behind using this approach is that 



dissipation along a given fluid element's trajectory occurs only at the shocks (unlike azimuthal 
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acceleration and deceleration due to pressure gradients that occurs also between the shocks) 
between the shocks the flow is adiabatic (in the absence of numerical dissipation). 



For weak shocks Larson ( 19901) gives the amount of energy irreversibly turned into heat at the 



shock per unit mass of the fluid crossing the shock as 



dE 
dm 



2 (M 2 - l) 3 



3(7 + 1) 



Mf 



where Landau & Lifshitz (1959) 



1 



7 + 1 AE 



(B3) 



(B4) 



4 E 

is the Mach number of the weak shock, AE/E is the density contrast across the shock, and s is the 
sound speed. For an isothermal shock (7 = 1) the expression M s = + (AE/E) is exact. Thus, 
for 7 = 1 we have 



dE_^( AE\" 



(B5) 



Now, in an annulus of width dzu the mass dm = 2m [Tt p — wY^dtdw crosses one of the 

2m shocks (the factor of 2 comes from accounting for both inward- and outward-propagating waves) 
per time d£, so that 



d w — 2m \ftjj — voYi- — , 
ow dm 



meaning that dF w ^ s /dzu = 2mwT,(dE/dm). Finally, 



/ dlY 1 



m 



:(M 2 



Ml 



m Es 2 



1 + 



AE\" 2 / AE\ 3 / din/ 



dlnzu 



(B6) 



(B7) 



If we define the "local disk mass" to be Yaw 2 then the accretion timescale t acc = Y>m 2 / M becomes 
(for AS/S < 1) 



_! 3 M 2 f AE 



rn w 



din/ 

dmcu' 



(B8) 



where M — s is the Mach number of the flow. 



If the disk were a M = const disk then M = 3na e ff s 2 E/f2, where a e fi is the "effective" 
a-viscosity due to trapped shocks. Combining this with equation (|B7[) one gets 



rn 



«eff 



din/ \ _1 (M 2 - 1) 



9tv V dlnzu 



(B9) 



in agreement with Larson (1990). 



